Astronomy & Astrophysics manuscript no. rmsynth3d 


©ESO 2011 


December 20, 2011 





Faraday synthesis 

The synergy of aperture and rotation measure synthesis 

M.R. Belffand T.A. EnBlirP 



o 

O 

Q 

00 



6 

c3 



> 
in 

<N 



Max Planck Institute for 
mrbell@mpa-garching . mpg . de 

Received ??? / Accepted ??? 



Astrophysics, Karl-Schwarzschild-Str. 1, 85741 Garching, Germany e-mail: 



Abstract 



We introduce a new technique for imaging the polarized radio sky using interferometric data. The new approach, which we call 
Faraday synthesis, combines aperture and rotation measure synthesis imaging and deconvolution into a single algorithm. This has 
several inherent advantages over the traditional two-step technique, including improved sky plane resolution, fidelity, and dynamic 
range. In addition, the direct visibility- to Faraday-space imaging approach is a more sound foundation on which to build more 
sophisticated deconvolution or inference algorithms. For testing purposes, we have implemented a basic Faraday synthesis imaging 
software package including a three-dimensional CLEAN deconvolution algorithm. We compare the results of this new technique to 
those of the traditional approach using mock data. We find many artifacts in the images made using the traditional approach that are 
not present in the Faraday synthesis results. In all, we achieve a higher spatial resolution, an improvement in dynamic range of about 
20%, and a more accurate reconstruction of low signal to noise source fluxes when using the Faraday synthesis technique. 

Key words. Polarization - Techniques: interferometric - Techniques: polarimetric - Methods: data analysis - Magnetic Fields 
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1. Introduction 

Rotation measure (RM) synthesis, introduced by 
iBrentjens & d e Bruyn (I2005I) . is a technique that makes 
use of the Faraday effect to improve the sensitivity of polari- 
metric observations by combining data over wide ranges in 
frequency. RM synthesis allows for the separation of polar- 
ized sources along the line of sight (LOS) by decomposing 
the observed polarized emission into parts originating from 
different Faraday depths (in the simplest case, Faraday depth is 
equivalent to RM), allowing one to generate a 3D representation 
of the polarized sky. While the Faraday depth axis cannot be 
mapped to a physical depth, the Faraday depth information 
can be of significant scientific value since the Faraday depth 
traces the projected strength and orientation of magnetic fields 
along the LOS. A more detailed introduction to RM synthesis is 
provided later. 

The RM synthesis technique has only recently become vi- 
able, owing to the availability of broadband receivers in the 
next generation of radio telescopes such as the Expanded 
Very Large Array (EVLA), the upgraded Westerbork Synthesis 
Radio Telescope (WSRT), and the pathfinder projects lead- 
ing to the Square Kilometer Array (SKA) such as the Low- 
Frequency Array (LOFAR). Recently, there have been many 
successful applications of RM synthesis, and interest is rapidly 
increasing as new radio telescopes are being commissioned. 
Applications have includ ed studies of the diffuse polarized emis- 
sion in the Perseus field (de Bruvn & Brentiens 2005; Brentiens 
I2011h and the Abell 2255 field ([Pizzo et al.l l2010h . analysis 
of the polarize d emission in ne arby galaxies in the WSRT- 
SINGS survey dHeald et al.ll2Q09h . and the detection of a shell 
of c ompressed magnetic fields surrounding a local HI bub- 
ble (IWolleben et al.l l201Qh . The RM synthesis technique will 
play a critical role in several upcoming polarization surveys, 



e^POSSUM (IGaensler et al]l2010h . GMIMS (IWolleben et al.l 
2009), and future surveys with LOFAR. RM synthesis i s very 
useful for studying magnetism. For instance, Be ll et al.l (l201ll) 
have shown that prominent asymmetric features in RM synthe- 
sis images known as Faraday caustics, which are related to LOS 
magnetic field reversals, can be used to study the structural and 
statistical properties of magnetic fields. 

In addition to the considerable interest in applications, there 
has been recent interest in improving RM synthesis imaging 
techniques. The RMCLEAN d econvolution algorithm was in- 
troduced by iHeald et all (I2009I) and was quickly adopted, ow- 
ing to its simplicit y and similarity to techniques used in aperture 
synthesis imaging. [Frick et all (1201 Ol) have proposed a wavelet- 
based RM synthesis technique. With this approach, one obtains a 
decomposition of the size scale of structures in addition to their 
Faraday depth location. There has also been gro wing interest in 
applying compressed sensing to RM synthesis dLi et al.l [20111: 
lAndrecut et al.ll201 ll) . Compressed sensing is a method of recon- 
structing signals that are sparse in some set of basis functions. 
If the signal is sufficiently sparse, it can be reconstructed us- 
ing fewer measurements than indicated by the Nyquist- Shannon 
sampling theorem. Since compressed sensing techniques make 
use of wavelet bases, but are implemented in a noise-aware fash- 
ion, they can be expected to be superior to pure wavelet based 
methods. 

These new imaging techniques have thus far focused on 
the problem of one-dimensional (ID) reconstruction. However, 
the product of RM synthesis is a function not only of Faraday 
depth, but also of position on the sky, making its reconstruction 
an inherently three-dimensional (3D) problem. In many cases, 
RM synthesis is performed on sky brightness images that have 
been produced from radio interferometric data. Observations 
performed with a radio interferometer sample the aperture plane 
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rather than the image plane, and this is done at many different 
frequencies. This data space is the 3D Fourier space represen- 
tation of the polarized sky brightness as a function of Faraday 
depth. Imaging algorithms should ideally make use of the en- 
tire data space to inform the reconstruction at each pixel, since 
the information about the sky brightness at each pixel is spread 
throughout Fourier space. However, the current approach is to 
perform the imaging in a piecewise fashion, first reconstructing 
the 2D sky plane images at each frequency before doing ID RM 
synthesis imaging along each LOS. Therefore, each step of the 
traditional imaging approach is done with a limited subset of 
the data, which will reduce the overall sensitivity and degrade 
fidelity in the final image. 

In this paper, we introduce a new technique for imaging the 
Faraday spectrum directly from radio interferometric data. We 
call this technique Faraday synthesis. Using this approach, one 
images the polarized emission as a function of sky position and 
Faraday depth from the visibility data itself, rather than using the 
traditional piecewise prescription. Faraday synthesis is a natural 
extension of the aperture synthesis plus RM synthesis techniques 
that provides improvements in image fidelity and sensitivity. 

We n o te tha t a similar technique was briefly discussed in 
Pen et al ] (120091) . although it was not considered in any detail, 
nor was it compared to the traditional approach to RM synthesis 
imaging. Furthermore, deconvolution was not considered. 

With the advent of RM synthesis the concept of rotation mea- 
sure, defined to be the amount that the observed polarization 
angle changes as a function of frequency, has become some- 
what outdated. With RM synthesis, one does not measure RMs, 
but instead reconstructs the polarized intensity as a function of 
Faraday depth. In the simplest case, where a single, discrete 
source of polarized emission is positioned behind a Faraday ro- 
tating medium, the RM is equal to the Faraday depth. In all 
other cases this is not true. In general, RM cannot be used as 
a proxy for Faraday depth, and the full distribution of polarized 
brightness as a function of Faraday depth is the most appropri- 
ate quantity to study. Therefore, we avoid use of the term RM to 
describe this new method, and instead call it Faraday synthesis. 
Throughout the remainder of this paper, RM synthesis will refer 
to the LOS imaging method developed by llBrentjens & de Bruynl 
( 2005). The traditional practical approach of first imaging indi- 
vidual frequencies using 2D aperture synthesis techniques and 
then reconstructing the LOS brightness distribution on a pixel- 
by-pixel basis will be referred to as 2+ ID imaging, in contrast to 
Faraday synthesis, which we will often refer to as 3D imaging. 

In Sec. [2] we briefly review the theories of aperture and RM 
synthesis imaging. In Sec. [3] we introduce the Faraday synthesis 
imaging technique. In Sec. [4] we describe the proof of concept 
software that we have implemented, and in Sec. [5] we compare 
test results obtained by imaging mock data using both the 3D 
and 2+ ID techniques. We conclude in Sec. [6] with a summary 
and discussion of our results. 



2. Synthesis imaging 

In this section we review the fundamentals of aperture and RM 
synthesis imaging before showing how they can be performed 
simultaneously in Faraday synthesis imaging. Reviews are in- 
cluded here for completeness and to highlight the assumptions 
that are typically made and the limitations that result. 



2.1. Aperture synthesis 

We can not possibly provide a complete review of the theory of 
aperture synthesis. We only wish to review those aspects that are 
most relevant to the current work. For a comprehen sive treat- 
ment, the reader is referred to lThompson et al.l (l200ll) . 

With an interferometer, one measures not the sky brightness 
directly, but rather a collection of discrete samplings of the aper- 
ture plane. These samples are complex quantities, typically re- 
ferred to as visibilities, denoted as V. The visibilities are the cor- 
related voltage output of pairs of antennas. For a narrow-band 
observation, they are related to the sky brightness distribution, /, 
by 

oo 

/dl dm 
1(1, m, A) 
Vl - P - m 2 

— oo 

e -2m[ul+vm+w( Vl-/ 2 -m 2 -l)] 

The coordinates (u, v, w) are spatial frequency coordinates, or 
the distance between pairs of antennas, measured in numbers of 
wavelengths. The coordinate u measures the distance in the car- 
dinal North-South direction, while v measures the distance in the 
East- West direction. The coordinate w points in the direction of 
the phase reference position on the sky. The (/, m) coordinates 
are direction cosines relative to the (u, v) coordinates. The wave- 
length at which the visibilities are measured is given by A. 

This relationship can be simplified to a two-dimensional 
(2D) Fourier transformation in two circumstances. The first is 
in the case of an East- West oriented array such as the WSRT. In 
this case, the telescopes move through a plane such that w = as 
the Earth rotates. The second case is when only a small patch of 
the sky is being imaged, such that w ( V 1 - P - m 2 - l) « 0. For 
the time being, we assume that we are looking at a small patch 
of the sky and henceforth neglect this w-term. 

A radio telescope is not equally sensitive to the entire sky. 
The sky brightness distribution is attenuated by the antenna 
power pattern, A, which is commonly referred to as the primary 
beam. Including this effect, the visibilities are related to the sky 
brightness distribution via the relationship 



Y (u,v,A) = jj dl dm A(l, m, A)I(l, m, A) ^" 2 * +vm) . (2) 

— oo 

In reality, as mentioned above, only discrete locations in the 
aperture plane are sampled. The measured visibilities, V, are re- 
lated to the true visibilities by 

V(u, v,A) = S (u, v, X)V'(u, v, A). (3) 
The sampling function S can be represented as 

S (u, v, A) = W(u, v, A) ^ 6 \u - j 

slv-^-\s(A-A i ) (4) 

where b = b x x+b y y is the distance between two antennas, known 
as the baseline length, and the unit vectors x and y point toward 
the North and East, respectively. The function W allows for the 
inclusion of weighting factors, e.g. by the inverse of the noise. 
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The i subscript is an index over the list of discrete values of b 
and A for which measurements have been made. 

To recover the sky brightness distribution from visibility 
data, one must invert Eqs. [2 and [3] Due to the sampling func- 
tion and the presence of noise, it is not possible to solve for / 
uniquely. The inverse Fourier transform of V does not give the 
sky brightness distribution alone, but rather 



m,A) = // du dv V(u, v, A) e 



2m(ul+vm) 



du dv S O, v, A)V\u, v, A) e 2nKul+vm) 
=B s Ul, m, A) * [A(/, m, A)I(l, m, A)] , (5) 

lm 

where L denotes convolution in the / and m plane. The image /d 
is commonly referred to as the dirty image. The dirty beam, B&y, 
is 



B s ky 



du dv S (u, v, A)e 



2m{ul+vm) 



(6) 



Due to the extended structure of the dirty beam, the sky im- 
age obtained from Eq. \5\ has a limited dynamic range and an 
unphysical appearance. Deconvolution is required to recover a 
reasonable approximation of / and to detect weaker features 
that are obscured by artifacts associated with the bright sources. 
By far, the most commonly used deconv olution a l gorith m is 
the CLEAN algor ithm, first intr oduced by Hogbom (1974) and 
later improved bv lClarkl (Il980h as well as lSchwabl (1 19841) . We 
describe the CLEAN algorithm later when discussing the 3D 
implementation that has been included in our proof-of-concept 
software. 

In the simplest case, when A is independent of baseline and 
time, the effect of the primary beam can be removed by sim- 
ply dividing by a known beam pattern. In doing so, the flux 
scale will be normalized across the sky and the noise level will 
increase as a function of distance from the pointing center. In 
general, however, A (and other so-called direction dependent ef- 
fects) can depend on both baseline and time, and in this case one 
mus t use something like the A-projection algorithm described 
bv lBhatnagar et alJ (120081) . 

2.2. RM synthesis 

Faraday rotation is a birefringence effect where the plane of 
polarization of a plane-polarized wave is rotated as it passes 
through a magneto-ionic medium. The right- and left-circularly 
polarized components of the plane wave propagate at different 
speeds through the medium causing a relative phase shift, and 
therefore a rotation of the polarization plane. The amount of ro- 
tation incurred is given by 



XU 2 ) = Xo + cpA\ 



(7) 



where x is the observed position angle at wavelength A, and xo 
is the position angle at the source of emission. The quantity 0, 
known as the Faraday depth, is given by 



z 

4>{z) = (0.81 rad/m 2 ) j |^ 



n e (z') 



cm - 



B z (z') 



(8) 



where B z is the LOS component of the magnetic field, and z is 
the distance along the LOS. The number density n e includes both 
thermal electrons and positrons, which are given as negative and 
positive values, respectivel y. Faraday d epth is distinct from RM, 
which we define following iBurnl (Il966l) to be 



(9) 



In the simplest case, when a single point source of polarized 
emission sits behind a Faraday rotating medium, the RM mea- 
sured for the source is equal to 0. In all other cases, these quanti- 
ties differ. In the general case of mixed Faraday rotating and syn- 
chrotron emitting media, the observed polarized intensity orig- 
inates from a range of Faraday depths, and the RM varies as 
a function of A 2 . The full polarized intensity as a function of 
Faraday depth, as obtained using RM synthesis, is required for 
an opportunity to study the intrinsic properties of the various 
sources along the LOS. 

The polarized intensity as a function of sky position and 
wavelength is a complex quantity that is usually defined in terms 
of the Stokes parameters to be 



P(l, m, A 2 ) = Q(l, m, A 2 ) + iU(l, m, A 2 ). 



(10) 



An important insight from Burn (1966) is that this is related to 
the polarized intensity as a function of Faraday depth, F, by 



P(l,m,A z ) 



oo 

/ 



d(p F(l, m, 



A 2 ) e 



,2i(f)A 2 



(ID 



We refer to F as the Faraday spectrum. The complex phase term 
reflects that the position angle of the polarized emission origi- 
nating at every location is rotated according to Eq.[7J 

Equation [TT] is similar to a Fourier transformation, and the 
ability to invert this relationship is desirable, but two things pre- 
vent this. First, it is not possible to sample P for all values of A 2 . 
One can only achieve a limited coverage of A 2 > 0, and of course 
cannot measure A 2 < at all. This pr oblem is addressed by 
the R M synthesis technique introduced by Brentjens & de Bruyn 
(120051) . where the inversion of Eq. [TT] is treated in much the 
same way as the inversion of Eq. [5j Second, the spectral depen- 
dence of F prevents one fro m inverting Eq. [TTJ As addressed by 
iBrentjens & de Bruynl (120051) . the inversion is possible assum- 
ing that the A 2 and dependent parts of F are separable, which 
means that the emission at all values along a LOS is produced 
with the same frequency spectrum (up to normalization), i.e. 



F(Z, m, 0, A 2 ) = /(/, m, 0>(/, m, A 2 ). 



(12) 



In general, of course, this assumption is not valid. Consider the 
case where two sources lie along the line of sight, each hav- 
ing different emission spectra. By factoring the Faraday spec- 
trum as above we introduc e an error into the image of /. 
IBrentjens & de Bruynl (120051) point out that such errors do not 
effect the Faraday depth of a source, and only have a relatively 
minor effect on the flux density. In their simulations, using a 
bandwidth that was 17% of the central frequency, an absolute 
error of 1 in the spectral index corresponds to an error of less 
than 5% in flux density. This error should be more pronounced 
with increased bandwidth, since the change in flux density over 
the frequency range becomes greater. 
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Assuming Eq.[l2] applies, the dirty Faraday spectrum, f D , is 
recovered by 

CO 

J s(l, m, A 1 ) 



=B^(l, m, (p) * /(/, m, (p) 



(13) 



where J is a convolution with respect to (p. The sampling func- 
tion, S A 2, is defined to be 

S A i(L m, A 2 ) = W#(U m, A 2 ) ^ S(A 2 - A 2 ), (14) 

i 

where A 2 are the discrete values of wavelength at which mea- 
surements have been made, and Wjz is a weighting term similar 
to that in Eq. |4] The so-called RM spread function, or the dirty 
beam in 0-space, B^, is given by 



£ (/,m,0)= dA z S A 2(l,m,A z )e 



,2i(pA 2 



(15) 



Equation [13] shows that, like with aperture synthesis, the 
product of RM synthesis is a convolution between the true 
brightness distribution (the Faraday spectrum in this case) and a 
dirty beam or point-spread function. Again, this dirty beam has 
structure that extends well beyond the source location. Therefore 
deconvolution is necessary in order to recover faint sources in the 
Faraday spectrum and obtain more physically realistic results. 

2.3. 2+1 D imaging of the Faraday spectrum 

RM synthesis is a ID imaging procedure that is applied to each 
LOS independently. With aperture synthesis imaging, one ob- 
tains the sky brightness distribution across the plane of the sky 
at a single frequency. Applying RM synthesis to such images 
requires that the variation of the sky brightness as a function 
of frequency must be determined at every location in the image 
plane. 

The sampling of the wv-plane, S , varies as a function of fre- 
quency. The images obtained from Eq.[5]at different wavelengths 
will be sensitive to different spatial frequencies and have dif- 
ferent resolutions. This complicates measurement of the spec- 
tral properties of the sky brightness because changes due to the 
sampling function can be confused with real variation of the sky 
brightness distribution. This problem is approximately overcome 
by applying different weighting and wv-plane tapering schemes 
to the data such that the wv-coverage is made to be roughly the 
same at each frequency. 

The usual practice of applying the RM synthesis technique 
to interferometric data involves several steps. We refer to this 
process as 2+ ID imaging: 

- Calibrate the visibility data. 

- Weight and taper the data such that the resolution of the im- 
ages is roughly constant with frequency. Some compensation 
is also needed if there is significant flux missing due to the 
usual gap in the center of the wv-plane. 

- Compile a series of deconvolved Stokes Q and U images at 
each frequency. Deconvolution is almost always performed 
using the CLEAN algorithm and the CLEAN model com- 
ponents are convolved with a so-called restoring beam, i.e. a 
Gaussian profile representing the resolution of the image (de- 
termined from the main peak of the dirty beam). The same 
restoring beam is used for all frequencies. 



- Stack the images. Pixel-by-pixel, the polarized intensity as 
a function of frequency is read from the maps, corrected for 
spectral variation, s (determined by measuring the spectral 
index from total intensity maps), and Fourier transformed. 

- Perform further 0- space deconvolution. 

The result is a 3D cube of data representing f(l, m, (p). 

There are two immediate problems with this procedure. First, 
the necessity to down weight data to approximately match resolu- 
tions between images at different frequencies can lead to signif- 
icant problems. This cannot be done perfectly and any variation 
in the polarized intensity arising from the differing resolutions 
produces a shift in the Faraday depth of the emission, thereby 
introducing systematic errors into the process. Second, Faraday 
synthesis is performed on maps that have already been processed 
using non-linear, ad-hoc deconvolution algorithms. Artifacts in- 
troduced into the images by said algorithms will be compounded 
during RM synthesis. 

We now introduce a new approach that is a natural extension 
of the typical synthesis imaging procedures described above, and 
does not suffer from the problems of the 2+ ID imaging tech- 
nique. 

3. Faraday synthesis 

In this section we show that it is possible to directly relate the 
Faraday spectrum to the visibilities of the linearly polarized 
emission. First, we decompose the Faraday spectrum into parts 
that relate directly to the intensity distributions of the two Stokes 
parameters Q and U(l, m, A 2 ). 

Equation [TT] can be rewritten as 



CO 

I 



Q + iU = d(p (F Q + iFu) e 



2iA 2 (p 



(16) 



where we have decomposed F into two complex valued terms, 
Fq and Fu, that relate directly to the Stokes Q and U brightness 
distributions, respectively. Therefore, 



Q = J dcpF Q e 2iA2 t 



(17) 



and the same relation holds for Stokes U, as well. We note that 
because Q and U are real, Fq and Fu are generally complex and 
Hermitian in (p, i.e. Fq(1, m, 0) = Fq(1, m, -cp). We note further 
that Fq and Fjj are not the local Stokes Q and U brightnesses 
at location (/, m, (p), but simply auxiliary variables useful for the 
formalism. If we assume that F can be factored as in Eq.[l2] then 



Q(l,m,A) =s(l,m,A) 



co 

J d<f>f Q (l, m,c 



2M 2 <f> 



=s(l,m,A 2 )q{l,m,A 2 ). (18) 

We will now work only with Stokes 2, but we note that the 
following expressions also hold for Stokes U and Fjj. Following 
Eq.[3] the measured Stokes Q visibility, Vq, is 



Vr 



oo 



dldmAQe 



-2m(ul+vm) 



=S JJ dldmAsqe- 2ni(uUvm) 



(19) 
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where S is defined as in Eq. H The relationship between the 
Stokes visibilities and the correlator output depends on the type 
of antenna feeds that are used. For linearly polarized feeds 

Vq =Vxx - Vyy 

Vu =V XY + Vyx, (20) 

where Vxx, Vyy, etc. represent the visibilities from cross- 
correlations of feeds having X or Y perpendicularly oriented 
dipoles. For circularly polarized antenna feeds, 

Vq =V rl + V LR 

Vu = - /(V RL - Vlr), (21) 

where Vrl and Vlr represent the visibilities from the cross- 
correlation between feeds of right and left, or left and right cir- 
cular polarization, respectively. 

We now define ail, m, 0) and cril, m, 0) to be the represen- 
tations of the primary beam and spectral dependence of the 
Faraday spectrum in Faraday space, respectively, i.e. 

oo 

A(Z, m, A 2 ) = J d</> a(l, m, ^)e 2iA ^ (22) 

— oo 

and 

oo 

s(l,m,A 2 ) = J d(Pa(l,m,(P)e 2a2<l) . (23) 

— oo 

Using the convolution theorem in between (p and A 2 , we combine 
the expressions above to find that the Faraday spectrum is related 
to the measured visibilities by 

oo 

Vq = S jjj dldmdcp (a*cr* /e ) e - 2 4'+™-^ (24) 

— oo 

where J denotes a ID convolution with respect to cp along each 
LOS. This expression can be inverted to give the dirty image for 
the Faraday spectrum 

oo 

{a * cr * /q j = JJJ du dv dA 2 y^ e 2m { ul+vm * ^) 

— oo 

=B * (a* cr * fo), (25) 

where 

oo 

B(l, m,<f>) = jjj dudvdX 2 S(u, v, ^ )e H ul+vm -^), (26) 

— oo 

and * is a full 3D convolution in /, m, and 0. We note that the 
3D and ID convolution operations do not commute. We can see 
that by inverting the visibilities, we do not directly recover the 
Faraday spectrum, but rather the Faraday spectrum convolved 
with B in 3D, and cr, and a in ID along each LOS. In order 
to recover /q, one must first perform a 3D deconvolution using 
the CLEAN algorithm, for example. After the 3D deconvolution, 
deconvolution of a and cr can be achieved by performing a ID 
inverse Fourier transform into /l 2 -space along each LOS, divid- 
ing by A and s, and then Fourier transforming back into 0-space. 

The beam pattern A is usually known to high precision and is 
often represented by an analytic function parameterized in /, m, 



and A 2 . In addition, a map of s will be required. This can be ob- 
tained by measuring the spectral variation of total intensity maps 
along each LOS. In many circumstances, it may be sufficient to 
simply assume that s is independent of / and m, since the errors 
introduced by using the wrong form for s are often quite small 
as previously discussed. 

Imaging software that implements the 3D inversion given by 
Eq. [25] would avoid the complications described for traditional, 
2+ ID imaging. We eliminate the need to match wv-coverage at 
all frequencies because the 3D dirty beam, B, is constructed from 
the full 3D sampling function. We also avoid the possibility of 
compounding errors through the process of deconvolving images 
that have already once been deconvolved. As a result, the fidelity 
of the images produced using Faraday synthesis should be im- 
proved over those made using the 2+ ID technique. In principle, 
the 3D approach will result in images that have higher dynamic 
range than those obtained using the 2+ ID approach because we 
are able to use all data across the full bandwidth during imaging 
and deconvolution. 

We have presented a simplified description of the Faraday 
spectrum measurement process above. This will already work 
quite well in many circumstances, notably for narrow-field ob- 
servations without significant direction dependent effects, but in 
general additional steps will be required. For instance, when the 
w-term in Eq .[T ] can no t be ignored, the w-projection algorithm of 
ICornwell et al.l (120051) has been shown to be very effective in re- 
ducing imaging errors. This algorithm makes use of the fact that 
the multiplication of the w-term and the sky brightness in the im- 
age plane is a convolution in visibility space. The w-dependent 
visibilities are projected onto the w = plane using the con- 
volution kernel, and then a 2D Fourier transform can be used 
to recover the sky brightness distribution. This algorithm can be 
applied to the case of Faraday synthesis without modification. It 
should also be possibl e to apply the A-projection algorithm of 
Bha tnagar et al.l (120081) . which corrects for direction-dependent 
beam effects in a manner similar to the w-projection algorithm. 

4. Proof of concept implementation 

To compare the 3D approach to polarization imaging with the 
traditional 2+ ID approach, we have implemented a proof of 
concept 3D imaging and deconvolution software package called 
f simager. For deconvolution, a 3D CLEAN algorithm has been 
implemented because it is by far the most common deconvo- 
lution method used in radio astronomical imaging, and with it 
we can make the most direct comparison between the two tech- 
niques. 

The CLEAN algorithm (Hogbom 1974) is a non-linear, it- 
erative deconvolution routine. The algorithm makes the implicit 
assumption that the sky is composed of point sources distributed 
throughout a mostly blank field. Over the last decades, it has 
been shown to work quite well even for fields that do not strictly 
meet this criterion. In brief, the procedure calls for iteratively 
building up a model of the sky by locating the peak of the im- 
age, adding a point source to the model at the location of the peak 
and with some fraction of its strength, and subtracting from the 
image the new model point convolved with the dirty beam. This 
is repeated until one can add no further flux to the sky model, 
i.e. when one is CLEANing the noise. A complete description of 
the 3D CLEAN algorithm that we have implemented is given in 
Appendix lAl 

Radio astronomical imaging relies heavily on Fourier trans- 
forms, and because the number of pixels in the 3D image of 
the Faraday spectrum is large, fast Fourier transforms (FFTs) 
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Figure 1. The sky-plane distribution of the 30 point sources in 
the model. The color scale indicates the Faraday depth, given in 
rad/m 2 , and the area of each circle is proportional to the log of 
the polarized flux. 



are required for computational feasibility. Visibility data is not 
collected on a regularly spaced grid in (w, v, /i 2 )-space. To use 
the FFT algorithm, the data must be interpolated onto regularly 
spaced grid points prior to processing. For this we employ a well 
known procedure known as gridding that is used extensively 
in aperture synthesis imaging as well as in medical imaging. 
Specifically, we ha ve implemented an algorithm described by 
Beat tv et al.1 (120051) . In brief, we convolve the data with a Kaiser- 
Bessel window function (KBWF) and sample the result on a reg- 
ular grid in (u, v, /l 2 )-space. After this procedure, one is able to 
perform a FFT with the same result achieved by using a discrete 
Fourier transformation, to within an arbitrarily small accuracy. 
The attenuation of the image plane caused by the convolution 
with the KBWF is corrected for by dividing the dirty image by 
the Fourier transform of the KBWF. There are a two parame- 
ters in the gridding procedure that affect precision at the cost of 
computational time. One parameter describes the extent of the 
KBWF in visibility space. The second parameter, the so-called 
oversampling ratio, is the factor by which the image plane should 
be enlarged in order to mitigate problems that occur at the edges 
of the image. With the parameters that we have chosen, for the 
KBWF to extend over 6 pixels in each of u, v, and A 2 , and an 
oversampling ratio of 1.5, the results are the same as one would 
obtain with a discrete Fourier transformation (DFT) to within 
one part in 10,000. In fact the dynamic range is only so limited 
near the edges of the image where the effects of the convolution 
with the KBWF are the most dramatic. Away from the edges of 
the image, the dynamic range is much higher. Overall, the dy- 
namic range can easily be improved by changing the gridding 
parameters, but this is done at the expense of increased process- 
ing and memory requirements. 

We have implemented fs imager in Python, thus allowing 
for rapid development and testing, but resulting in overall poor 
performance. To improve the situation, we have optimized some 
sub-functions (particularly the gridding routines) using Cython 
and in the end the code performs admirably. We can load a 1GB 
data set, grid, image, and CLEAN a 16 megapixel image using 
500 iterations in about 15 minutes on our 2.4 GHz Core i5 devel- 
opment machine with 8 GB of RAM. The most major limitation 
of the current version of the software is that all of the data re- 



sides in memory and this limits the size of the images that can 
be produced to the amount of memory that is available on the 
machine. We have run all tests on a computer having 64 GB of 
memory, and are limited to producing image cubes that are 400 
megapixels in size (about 750 pixels per side). This may be suf- 
ficient for imaging small fields of view, but when imaging data 
from a wide-field, high resolution instrument (e.g. LOFAR), this 
is inadequate. 



5. Tests 

To compare the 3D and 2+ ID approaches, we have produced a 
mock observation of a set of polarized point sources. In this sec- 
tion, we first describe the mock data, and then we describe the 
method that we use to image the data using the 2+ ID method. 
Finally, we compare the Faraday spectra produced by each tech- 
nique. 



5.1. Mock data 

We create a model Faraday spectrum consisting of a number of 
polarized point sources. We choose this simple scenario, though 
it may not be the most physically motivated one, because it is 
what the CLEAN algorithm was designed for. In these tests we 
simply want to study the differences between the two imaging 
frameworks, not the merits of the specific imaging algorithm 
being used. By choosing a model for which CLEAN is ideally 
suited, we can concentrate on the fundamental differences be- 
tween the 3D and 2+ ID techniques without being concerned 
about artifacts introduced by using an inappropriate deconvo- 
lution algorithm. 

No spectral variation is included in the model Faraday spec- 
trum, i.e. F is only a function of /, m, and 0. Also, we assume 
that our field of view is small relative to the size of the primary 
beam, and therefore that there is no variation in beam strength 
throughout the image plane. 

Our model includes 30 polarized point sources, randomly 
distributed throughout the volume. The distribution of model 
sources is shown in Fig.[T] The color indicates the Faraday depth, 
where the scale shown on the right is in rad/m 2 . The size of the 
points in the plot indicates the relative linearly polarized flux of 
the source. The fluxes range from 0.06 Jy to 64 Jy, thus allowing 
for a comparison across a wide range of signal to noise ratios. 

We "observe" the model by taking its Fourier transfor- 
mation and then filling the data column of a custom made 
MeasurementSet file. The MeasurementSet is created using the 
makems script that is part of the LOFAR software package. This 
script creates a MeasurementSet with user specified time and 
frequency coverages, and computes wv-coordinates using the an- 
tenna table from a pre-existing MeasurementSet. We use the an- 
tenna table from a VLA A-array observation. The frequency cov- 
erage is specified such that our mock data spans the range from 
1-4 GHz, with 64 frequency channels distributed throughout the 
range. The total time covered by the observation is 1 hour, with 
a 60 second step size between measurements to reduce the file 
size. 

Gaussian white noise, having a standard deviation of 10 Jy, is 
added to real and imaginary parts of the stokes Q and U visibil- 
ities separately. The MeasurementSet file is then either read di- 
rectly into the 3D software for gridding, imaging, and deconvo- 
lution, or into the CASA software package for traditional imag- 
ing. 
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Figure 2. The 3D dirty beam, a) A 2D slice through the cube in the sky plane at = rad/m 2 . The magnitude of the complex valued 
beam is shown. The bottom contour begins at 0.01 and there is a step of 2 between levels. Negative contours are marked with dashed 
lines, b) Another view of the beam, now at = 50 rad/m 2 . c) A ID profile through the center of the sky plane at (0", 0"). d) A ID 
profile along the LOS at (-3.5",0.2"), marked with an "x" in Figs, a and b. The magnitude of the profiles are shown as a solid line, 
while dashed and dotted lines indicate the real and imaginary parts, respectively. 



5.2. 2+1 D imaging procedure 

The 2+ ID imaging procedure is conducted by first loading the 
mock observation data into CASA where it is then written to 
UVFITS format using the task exportuvfits. This data is loaded 
into AIPS and each channel is independently deconvolved using 
the IMAGR task. 

In an attempt to match the wv-coverage at the different fre- 
quencies, weighting and tapering schemes are employed such 
that the output image resolutions are roughly equal. At 1 GHz, 
the maximum wv-spacing is approximately 125 kA and the main 
peak of the synthesized beam has a FWHM of approximately 
1.7". At 4 GHz, the maximum wv-spacing is 490 kA and the 
FWHM of the main peak of the synthesized beam is approxi- 
mately 0.6". Without tapering, we will observe large variations 
in the sources as a function of frequency simply due to the dra- 
matic change in resolution. 

All but the lowest frequency observations are tapered in the 
wv-plane with a Gaussian profile having a FWHM of 150 kA. 
Robust weighting is used, and the weighting parameter is cho- 
sen such that FWHM of the main peak of the synthesized beam is 
approximately 1.2". After deconvolution in the image plane, all 
images are restored using a Gaussian profile that has a FWHM of 
1.2". We have tried a few other weighting and tapering schemes 
so to achieve even lower resolutions, down to a FWHM of 1.7". 



The different schemes made no significant difference in the re- 
sults apart from the resolution of the final image cube. 

RM synthesis is performed along each LOS using our own 
RM synthesis software that is currently in use on the LOFAR 
compute cluster for commissioning. T his software imple ments 
the RMCLEAN algorithm described bv lHeald et aD (120091) . 

5.3. Results 

The result of either the 2+ ID or 3D imaging procedures is a 3D 
reconstruction of of F(l, m, 0). In Fig. [3] we show a side-by-side 
comparison of the reconstructed Faraday spectra for a few se- 
lected Faraday depths. In each row, the left panel in each image 
shows the model image, the middle panel shows the 2+ ID re- 
construction and the Faraday synthesis reconstruction is shown 
on the right. 

The most obvious difference between the two images is the 
resolutions. In the 3D reconstruction, we use a natural weight- 
ing scheme. The FWHM of the main peak of the 3D dirty beam 
is roughly 0.8", 0.8", and 40 rad/m 2 in R.A., Dec, and 0, re- 
spectively. Selected image plane slices, and LOS profiles of the 
Faraday synthesis derived 3D dirty beam are shown in Fig. 13 For 
the 2+ ID imaging, we have in some sense chosen the resolution 
in the sky plane to be 1.2" by selecting a particular weighting 
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Figure 3. A comparison between the input model and the two reconstructions at several cp values. In each frame, the greyscale varies 
linearly from 0-50 mJy/beam, and the contours begin at 50 mJy/beam with a factor of 2 between each level. Left: The model image 
convolved with an (0.8"x0.8"x40 rad/m 2 ) Gaussian. Middle: The 2+ ID reconstruction. Right: The fsimager reconstruction. Top 
to bottom:^ = -205, -160, 200, 230, and 395 rad/m 2 . 
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Figure 4. A few LOS profiles through the reconstructed Faraday 
spectra. Note that the brightness scale is logarithmic. In each, the 
solid and dashed lines indicate the 3D and 2+ ID reconstructions, 
respectively. The profiles are taken along the following image 
plane coordinates: a) (7.2", 3.2"), b) (5.3", 5.3"), c) (-6.5",-7.4"), 
d) (-7.9",-4.3"). 




io u 

Model fluxes (Jy) 

Figure 5. Reconstructed source fluxes compared to those of the 
input model. Dark circles denote the fsimager reconstruction 
while light crosses denote the 2+ ID reconstruction. Arrows in- 
dicate upper levels of reconstructed fluxes that were not located 
by our search algorithm. 



scheme, but our choice is limited by the resolution of the low- 
est frequency part of the data. We are therefore able to achieve a 
higher resolution with the Faraday synthesis technique, without 
the need for tapering the long baseline data. 

The LOS profile through the center of the 3D dirty beam, 
shown in Fig. Efc, is the same as the usual RM synthesis dirty 
beam given in Eq. Q3J However, the off-center profile, shown in 
Fig. [2fc, is not simply a scaled version of the usual RM synthesis 
dirty beam. The structure of the two profiles is quite different. 
This is because the sidelobes of B s ^y change as a function of 
frequency, which leads to structure in 0-space. 

We find that the noise level in the 2+ ID image is higher than 
that of the 3D image. We measure the noise levels in the two 
reconstructions by computing the root mean square (RMS) pixel 
values in several regions of the image cubes where no sources 
are located. In the 3D image cube the noise is 6.06 mJy/beam. 
The RMS is 20% higher in the 2+ ID image cube, measuring 
7.33 mJy/beam. 

Both imaging techniques are able to successfully reconstruct 
the stronger sources in the model image quite well. The sources 
above 1 Jy have all been located correctly, and the recovered 
fluxes are within a few percent of the model fluxes in each case. 
In the reconstruction at = -205 rad/m 2 , shown in the first 
row of Fig. [3j both the 2+ ID and 3D reconstructions roughly 
match the input model. Even the weak sources are detected. This 
represents one of the better slices of the 2+ ID reconstruction, 
and yet here we can already begin to see problems. There are 
some hints of artifacts in the traditionally made image that are 
not present in the Faraday synthesis result. The artifacts appear 
prominently in most other frames of the reconstruction. 

Overall, the blank sky in the Faraday synthesis reconstruc- 
tion is noise-like throughout the image volume, while in con- 
trast, we find many artifacts present in the 2+ ID reconstruction. 
The artifacts often appear as circular outlines of sources at dif- 
ferent Faraday depths, like the feature at ~(7", 3") in the image 
at cp = -160 rad/m 2 . Such features are present in each of our 
example image slices above. These circular features are quite 
obviously artifacts, but in some cases, for example at ~(-6", -7") 
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in the image at cp - 200 rad/m 2 , false features are quite strong 
and not as obviously artificial. One can see a weaker analog of 
this particular artifact in the 3D reconstruction as well, but most 
other artifacts in the 2+ ID images have no counterpart in the 
3D reconstruction. These erroneous features act to confuse the 
detection of weaker sources in the field. 

There are artifacts that appear along the LOS profiles of the 
2+ ID reconstruction as well. In Fig. [H we show a few LOS 
profiles through each image cube. Figure shows a profile 
through a relatively strong and isolated source at R.A.=7.2", and 
Dec. =3. 2" (c.f. Fig.Q]). This source is well reconstructed by both 
methods. In the other example profiles, problems with the 2+ ID 
reconstruction become apparent. The profile shown in Fig. 01) is 
again through a rather isolated but weaker source, at R.A.=5.3", 
Dec. =5. 3". Here we find that the actual source is not prominent 
compared with the noise in the 2+ ID profile, but is easily de- 
tected in the 3D reconstruction. 

The other profiles are through more crowded regions of the 
image cube, where contamination from nearby sources becomes 
problematic. The profile in Fig.|4t, along R.A.=-6.5" and Dec.=- 
7.4", passes through another relatively strong source. While the 
source itself is well reconstructed, three other features due to 
nearby sources appear around =-150, 200, and 350 rad/m 2 in 
the 2+ ID reconstruction. These are not present in the 3D case. 
Figure |4ji shows a profile through a relatively weak source at 
R.A.=-7.9", and Dec. =-4. 3". The source is easily distinguish- 
able in the 3D profile, but in the 2+ ID profile it is completely 
overshadowed by the artifact attributed to the source at R.A.=- 
6.7", and Dec. =-5. 7", almost 2" away. 

These problems occur in part due to the lower resolution of 
the 2+ ID reconstruction. They are also due to the fact that each 
channel in the data is imaged and deconvolved (in 2D) sepa- 
rately, and therefore with a more limited sensitivity. As a result, 
residual artifacts remain in the individual images, and these be- 
come apparent with the increased sensitivity provided by com- 
bining data over the full frequency range. These residuals are 
then processed during RM synthesis, and can become quite prob- 
lematic as we can clearly see in our examples. 

In each reconstruction, even the weakest model source 
should be present at the 10<x level in the 3D image and at 8cr in 
the 2+ ID case. To locate point sources in each image, we have 
a simple routine for locating local maxima that scans through 
every pixel and checks whether it is larger than all neighboring 
pixels. If so, the location and brightness is recorded. We search 
through each image cube in this way to find all local maxima 
above the 50 mJy/beam level. Ideally we should expect to find 
30 points corresponding to the model sources. In the 3D recon- 
struction, we find 32 such locations. All input model sources are 
located along with two false detections. In the 2+ ID case we lo- 
cate 147 sources, with 5 input model sources not detected. A few 
of the missing sources are not detected because they are simply 
not resolved from a nearby stronger source. The others may be 
present below the 50 mJy/beam cutoff, but at this point they are 
completely indistinguishable from the artificial sources. 

The artificial sources in the 2+ ID reconstruction are not 
only much more numerous, but also brighter. The brightest ar- 
tifact in the 2+ ID image cube is 0.113 Jy/beam, compared to 
0.062 Jy/beam in the 3D image cube. 

Figure \5\ compares the fluxes of the sources found in the two 
reconstructions to those of the input model. These have been cor- 
rected for the Ricean bias effect described in lWardle & Kjonbergl 
(119741) . The 3D reconstruction agrees remarkably well with the 
model fluxes across the full range of source strengths. This sky 
model is ideally suited for the CLEAN algorithm, so this is as 



should be expected. The stronger sources are also recovered 
nicely in the 2+ ID reconstruction, but the strengths of weaker 
sources are systematically lower than the input model fluxes. 

6. Discussion and conclusions 

We have described a new approach to imaging linearly polar- 
ized visibility data that we call Faraday synthesis. With this ap- 
proach, one directly reconstructs the Faraday spectrum, or the 
polarized intensity as a function of Faraday depth, from the vis- 
ibility data. This takes the place of the usual approach of first 
performing aperture synthesis imaging on the visibility data at 
each frequency, then performing RM synthesis along each line 
of sight independently. 

These two approaches would be equivalent if deconvolution 
were not required. With Faraday synthesis, the deconvolution 
is done in one step using the entirety of the broad-band data. 
In contrast, the traditional approach requires deconvolving the 
images individually at each frequency. The sensitivity in the 
narrow-band images is significantly limited, and residual arti- 
facts remain in these images that limit the dynamic range and 
image fidelity achieved during RM synthesis. Moreover, another 
deconvolution procedure is necessary to account for the point 
spread function due to the incomplete sampling of the wave- 
length space. Artifacts introduced by the first deconvolution al- 
gorithm will be compounded during this procedure, further re- 
ducing image fidelity. 

Indeed, we found in our proof-of-concept testing that arti- 
facts were significantly higher in the traditional imaging method 
than with Faraday synthesis. The noise was roughly 20% lower 
when using the Faraday synthesis technique, and the strongest 
artifact was about half as bright. 

We found that one is able to achieve a better resolution in the 
final image using the new approach. The main lobe of the 3D 
dirty beam was 30% smaller in the sky plane than that of the tra- 
ditional method. With the 2+ ID method, stacking of the images 
at each frequency requires tapering the visibility data such that 
the higher frequency images have the same resolution as those 
at the lower frequencies. We find that the inaccuracies inherent 
in the stacking process are a significant source of artifacts in the 
traditional RM synthesis technique. Furthermore, this procedure 
requires severely down-weighting a large portion of the avail- 
able data, again limiting sensitivity. With the Faraday synthesis 
approach, no such down- weighting is required. 

The Faraday synthesis approach of working directly between 
visibility and Faraday space is a much better foundation on 
which to build new image reconstruction algorithms because 
with it one is able to accurately describe the instrument response. 
The effects of the intermediate, non-linear deconvolution proce- 
dure can not be easily understood or modeled. Many signal in- 
ference techniques, l ike those derived usin g Info rmation Field 
Theory (EnBlin et al. 2009; EnBlin & Weig 2010), depend on a 
complete description of the instrument response and would need 
to be built on the framework of Faraday synthesis. 

While the CLEAN algorithm has worked quite well in radio 
astronomy for decades, the implicit assumption of a sky sparsely 
populated by point sources is not well suited for the kinds of 
diffuse polarized signals that one finds, for example, in the po- 
larized emission from the Milky Way. Novel image reconstruc- 
tion algorithms built using more appropriate constriants or as- 
sumptions are likely justified. Such algorithms could make use 
of statistical correlations inferred from the d ata, similar to the 
extend ed critical filter algorithm developed by lOppermann et ail 
(1201 ll) . 
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after the model is subtracted from the visibilities are called the 
residual visibilities, Vr. 

Before starting the procedure, we extract a patch from the 
dirty beam, # pa tch, for use during the minor cycle. We record the 
value ft = max(|# - # pa tchl), where | • | indicates that we take the 
magnitude of the complex valued map. 

The major cycle includes: 

1 . Invert Vr to create a residual image of the complex Faraday 
spectrum, Fr. Find Fn m = max(\F^\). 

2. Start the minor cycle using Fr and populate a minor cycle 
sky model, M^ nov . The minor cycle is described below. 

3. Upon completion of the minor cycle, inverse Fourier trans- 
form the minor cycle sky model into visibility space, Vm. 

4. Subtract the model from the residual visibility data, Vr = 
Vr-SVm. 

5. Add Mminor to the complete sky model, M ma j or . 

6. Repeat steps 1-6 until some user defined number of iterations 
has been done, or F\ im is below a user defined cutoff. 

7. Invert Vr to produce a final residual image. Add to this 
^major convolved with a Gaussian restoring beam. 

The minor cycle proceeds as follows: 

1. Find (/ m , m m , m ) = argmax(\F R \). 

2. Add a point source to Vm at (/ m , m m , cp m ) with a flux Fm = 
gFR(l m , m m , m ), where g is a user defined gain factor (be- 
tween zero and one). 

3. Subtract ^M^patch, centered on (/ m , m m , m ), from F R . 

4. Continue as long as \Fm\ > f3Fn m (l + ^) where N is the total 
number of minor cycle iterations that have been completed. 

We a dopt the minor cycle stop condition suggested in Clark 
(1980). This condition reflects the fact that during the minor cy- 
cle we only subtract a small patch of the dirty beam from the im- 
age and therefore any effects outside of this patch will not be re- 
moved from the image. During the minor cycle, we CLEAN only 
as deeply as the maximum contribution of the brightest source to 
the residual image that is not removed by subtracting the beam 
patch. Actually, the (1 + ^) term makes this condition even more 
strict early on in the procedure. 



Appendix A: 3D RMCLEAN 

In our proof of concept software we have impl emented a va riant 
of the CLEAN routine that was introduced bv lCla^([T98Qk . In 
this variant the model is populated by searching for peaks in 
image space, but the model subtraction is performed in wv-space. 

We first produce the dirty images of Fqv and Fjjq . These are 
combined into a single complex image according to 

%(F D ) =%(Fqp) - 3(Fud) 

3(F d ) =3l(F ro ) + 3(F0d), (A.l) 

where % and 3 are operators that select the real and imaginary 
parts of the image, respectively. We proceed with the CLEAN 
procedure using the complex valued Fd, and the associated vis- 
ibilities. 

The algorithm is performed in two parts, the so-called major 
and minor cycles. In the minor cycle, new model sources are lo- 
cated in image space. The subtraction of the model from the vis- 
ibility data is performed during the major cycle. What remains 
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